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*^ We present a model of spheres moving in a high- dimensional compact space. 

Ch We relate it to a mixed matrix model with a 0(N) invariant model plus 

o 

a P(N) invariant perturbation. We then study the low pressure regime by 
performing a diagrammatic expansion of this matrix model. Finally, we show 
the results from numerical simulations and we present evidence for a glassy 
regime at high pressures. 
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Glassy systems can be characterized as out-of-equilibrium systems in which the relaxation 
is slower than the experimental time-scale [1]. The slowness of the dynamics is due to the 
complexity of the energy (or free-energy) landcape, which may contain a combination of 
barriers, bottlenecks and even flat directions. In the case of spm-glasses, the complexity of 
the landscape is induced by the frustration in the spin configurations, which in turn is a 
consequence of the (quenched) disorder in the couplings. In structural glasses the ruggedness 
of the landcape does not necessarily come from a quenched disorder, but is a consequence 
of the geometrical constraints for the motion; the disorder is self-induced. 

Despite this difference in the origin of disorder, both spin and structural glasses share 
some common features such as a long-term memory going back to the time when the system 
was quenched to low temperature. These 'aging' effects can be studied analytically in simple 
mean-field spin-glass models [2] and they are qualitatively very similar to the aging effects 
found in true experimental spin-glasses. 

Recently there has been a considerable success in constructing the 'structural' counter- 
part of such simple spin-glass models, i.e. mean-field systems without quenched disorder 
having a dynamic phase transition into a glassy phase. These models include spin models 
[4-7], particles on a hypercubic cell [8] and a field theory with quasi-random interactions [9]. 
In all but the latter model, the glassy behaviour is a consequence of the discreteness of the 
variables: a continuous version of them behaves like a 'liquid' at all temperatures. 

The purpose of this paper is to present a model of hard spheres [3] in a high dimensional 
(off-lattice) compact space. The model is related to (but different from) ordinary spheres in 
a flat high-dimensional space [10]: as we shall show below it has a more complicated low 
pressure virial expansion, but it is easier to simulate for high dimensions. We obtain analyt- 
ical results for the low pressure regime, and numerical evidence for a glassy transition. We 
simulate the constant-pressure dynamics in the high-pressure phase, and study the ergodic 
properties in that case. 
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I. THE MODEL 



We consider a set of N particles of positions S a (a = l,...,iV) constrained to be on a 
D- dimensional sphere: 

^XK a ) 2 = l, Va (1) 

8 = 1 



The 'distance' between two particles on the sphere is: 

q ab _ 1 

and the hard sphere condition reads: 



Y^S?S! = cos 6 ab , Va<6 



-Q < q ab < Q , Va < b (3) 

where Q is the 'size' of the particles. We choose Q of order one so that the angle between 
particles 9 ab on the sphere is close to 90° (cfr. Eqs.(l),(2)). This yields a non-trivial be- 
haviour in the high-dimensional limit with the particle number scaling with the dimension 
as N = pD and the 'density' p of order one. For every particle there is another in the 
antipodes of the sphere. Equivalently, we can say that we identify opposite points on the 
sphere. In that case, the manifold on which the particles move is no more the D-dimensional 
sphere but the projective plane RP(D). 

An alternative normalization is to impose the minimal angle between any two particles 
Gmin to be fixed and different from 90° in the limit N — > oo. In that case, Q is of order \[~D 
and N has to grow exponentially with D to have a non-trivial behaviour [3]. Spheres in a 
flat high-dimensional space are recovered in the limit Q m in 0. We do not consider here 
this alternative normalization. 

In order to study a constant-pressure situation, one can consider a variable radius of the 
sphere on which the particles move, and a coupling term between the radius of that sphere 
and the pressure. We shall instead follow the equivalent procedure of keeping the radius of 
the spherical space fixed, and considering the particle size Q as a dynamic variable. The 
partition function reads 



r r N 

Z[p,P] = JdQ exp[-N 2 PQ] J ]JdS a 6 (S 2 a -D)l[6 (qV15 - \S a .S b \) . (4) 

a=l a<6 



The coupling with the pressure P is provided by the factor exp[— N 2 PQ] in the partition 
function. Note that any other function of Q just amounts to a (nonlinear) rescaling of the 
pressure. 

In the Section III we shall study the response of the density of the system to a space- 
dependent field. In an ordinary flat space one applies a field with a strength that varies 
spatially as hit) cos(kx) and then measures the corresponding space Fourier component of 
the density. Here we shall apply a field that pushes the particles towards the direction 



Ht) £ 4 (5) 

ia 

and we shall then measure the corresponding density response 

m(t) = < E *'«(*) ) • (6) 

Finally, let us remark that one can define a shear viscosity for this model by subjecting 
the system to a tangential force that varies with the 'latitude' from the equator of the 
spherical space, and measuring the corresponding average velocities. We do not do so here. 



II. ANALYTICAL CALCULATIONS 

In the thermodynamical limit, the free-energy depends only on the pressure P and on 
the density p = N/D. Introducing the entropy 




of all positions corresponding to a given set of overlaps q ab (Eq.(2)), Z[p, P] may be rewritten 
as the partition function of the following matrix model 

Z[p,P] = JdQ exp[-N 2 PQ] f n [dq ab e{Q - \q ab \)\ e s ^ (8) 
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The origin of the difficulty arising in the computation of the partition function is now 
clearer. Whereas any rotation of the matrix q ab in the N X N- dimensional space of the 
particle-overlaps leaves the entropy S(q ab ) unchanged, the hard sphere condition is only 
invariant under the group of permutations P(N) of the particles. As a consequence, the 
action in Eq.(8) does not depend only upon the eigenvalues of q ab and the angular variables 
cannot be integrated out to obtain a solvable model in the large N limit [14]. We shall now 
see that this difficulty does not prevent us from estimating the full diagrammatic expansion 
of the free-energy which is expected to be valid in the low P regime. 



A. The diagrammatic expansion 

Let us consider now a generic matrix model M ab whose action includes two different 
terms. The first one, S\(M) } is O(N) invariant. The second one, Y^ a b ^(Ma?,), breaks the 
rotational symmetry but is invariant under a permutation of P(N). Both parts S\ and S 2 
may be expanded in integer series of their argument. Assuming that the mean value of M ab 
is equal to zero, the first terms are quadratic in M ab while the following ones contain all 
kinds of even vertices compatible with the O(N) and P(N) symmetries : 

Si(M) ^-E -W: • y E M ab M bc M cd M da + 0(M 6 ) (9) 

a,b a,b,c,d 

S2(M) = |E M l + 92 £ Ml + 0(M 6 ) (10) 

a, 6 a, 6 

where the normalization factor in Eq. (9) ensures that the limit N — > 00 is well defined. For 
simplicity we shall consider only quartic interaction terms; the extension of what follows to 
the general case is immediate. We therefore end up with the computation of 

V jU'M ^vl ^J»: J E M ab M bc M cd M da -g 2 J2 -W ' ] (11) 

a,b \ a,b a,b,c,d a,b J 

where r = + r 2 . When g 2 = 0, Y may be expanded in powers of §\ and only planar 
diagrams survive in the large N limit [14]. If we now consider the diagrams coming from 
the g 2 vertex, the self-energy will include some new contributions which may be separated 
in two groups : 
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• some diagrams are topologically equivalent to the ones arising in the g 2 = expansion 
but the propagators between the §\ vertices are now "dressed" with the g 2 interaction : 
such diagrams may be simply obtained by the addition of a self-energy term r 2 + 
<J 2 J2 a b M 2 b m where the value of o 2 will have to be fixed in a self-consistent 
way. 

• the remaining "mixed" diagrams include interacting §\ and g 2 vertices : the simplest 
representative of this class is 

^ J2 (M ah M hc M cd M da Mt } ) (12) 

where (•) denotes the Gaussian bare measure over the matrix elements. Once Wick's 
contractions have been done in the inner loops of such diagrams, the number of free 
indices running inside the loops is strictly lower than the power of the j-j factor. 
Following the lines of [5], it seems reasonable to conjecture that this result still holds 
for all "mixed" diagrams. Therefore, these diagram should not contribute in the large 
N limit. Although we have not been able to prove this conjecture, we believe it is 
exact. 

The S 2 term turns out to renormalize the self-energy of the 0(N) invariant model by 
an additional factor o 2 . As a consequence the average squared value (M 2 ) equals the full 
propagator of the Si theory, (M 2 )i(r + <t 2 ) with the bare mass r + o 2 . 

If we now consider S\(M) as a perturbation of the P(N) invariant action S 2} the above 
reasoning will still hold and the full diagrammatic expansion is identical to the one of the 
S 2 model with a new contribution <Ti to the self-energy. Thus, the mean squared value of 
the matrix element equals (M 2 ) 2 (r + 01) where (-) 2 denotes the average value with the S 2 
action and the bare mass r + a\. It is furthermore equal to the inverse of the renormalized 
mass. The self-consistency equations that <Ti and o 2 fulfill are therefore 

(M 2 Mr + a 2 ) = (M 2 ) 2 (r + a,) = -— L (13) 

r + (Ji + cr 2 
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B. Resummation and introduction of the disorder 



We shall now see that the above diagrammatic expansion may be obtained in a compact 
form by introducing the new partition function (on the same lines as in Ref. [5]) 

Y[U] = J n dM ab exp M lb " 77 E M ab M bc M cd M da - g 2 J2 (e U ab , cd M cd \ 



\ 



a,b,c,d a,0 \ c,t 

(14) 



where U is a N 2 X iV 2 orthogonal matrix. The important point is that the annealed par- 



tition function ^[f], averaged over all orthogonal matrices U, has the same diagrammatic 
expansion as the original Y (11) which corresponds to the particular choice U abjCd = S ac 6 bd . 
Let us define T ab = Y^ c d U a b,cdM cd and its Lagrange multiplier T ab to rewrite 

Y[U] = jH<l.\l : <IT : <IT : vx\J ^ ^ J2 M ab M bc M cd M da 

a,b \ a,b a,b,c,d 

- 92j2 T ab + J2Ta b T ab - J2 T ab U ab , cd M cd \ (15) 

a, 6 a, 6 a,b,c,d J 

If we now average Y[U] over all possible choices of the orthogonal matrix £/, we obtain [5] 
Y\U] = jH<l.\l : <IT : <IT : vx\J ^ ^ J2 M ab M bc M cd M da 

a,b \ a,b a,b,c,d 

- 92 E T A ab + £ f ab T ab + g(j: T 2 ab £ K\ ) (16) 

a, 6 a, 6 y a, 6 c,d J J 

where 

G(t) = log J 1 dx{\- x 2 ) N l 2 e lNx ^ (17) 
Introducing the order parameters ->W Yla,b Tab an d ^ = ~Ap Yla,b an d their respective 



Lagrange multipliers r and A, we see from the above expression (16) for Y[U] that the 
integral over the T variables is purely Gaussian and may be analytically performed. The M 
and T variables, whose vertices are respectively 0(N) and P(N) invariant, are thus coupled 
only through an additive renormalization of their masses by A and r respectively. Therefore, 
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the introduction of the disorder and the annealed average we have carried out has permitted 
us to find another partition function which must be equal to the whole resummation of the 
diagrammatic expansion of Y exposed in the previous paragraph, provided that Eq.(13) is 
fulfilled, that is 

<M 2 )(A) = <T 2 )(r) = -±- (18) 

T + A 

It is possible that the disordered model also yields information about the high pressure 
phase. 

C. The hard spheres free— energy 

Let us now apply the previous results to our original problem of interacting hard spheres. 
Starting from Eq.(8), we introduce the orthogonal matrix U 

N 



Z[U, p,P] = j dQ exp[-N 2 PQ] j f[ dSJ (S 2 a -D)Y[el QV3 

a=l a<b \ 

and average over all possible choices of such matrices to find 



y ] Uab,cdS c -S a 
c<d 



(19) 



Z[U, p,Q} = j dQd\dre- N2F{x ' T ' Q ' p '^ (20) 

where 

F(X, r, Q, p, P) = XF P (X, p) + T -F s (t, Q) - i log(r + A) + QP (21) 

and the two free-energies corresponding to the 0(N X N) and P(N) interactions respectively 
equal 



-^log / l[dS a S(S 2 -D) exp f-^E(^V) 

a=l \ a^b J 



F p (X,p) = ^m^-— log I l[dS a S(S i a -D)e X p\-—Y:(S a .S b y\ (22) 
at fixed density p = and 

F s (t, Q) = -- log [ Q exp {--x 1 ) (23) 

T J-Q V27T V 2 

The saddle-point equations with respect to r and A read 



(^) p = <*>). = ^ CM) 

where the measures (.) p and (.) s correspond to the free-energies (22) and (23) respectively. 
Defining explicitly the latter, 

exp [ —t:X 2 ) 

V(x) = — V 2 > (25) 

Lq dy exp \ -\y 2 ) 

the saddle-point over Q reads 

V(Q) = P (26) 

Interpreting V(x) as the probability distribution of the overlap between two spheres, identity 
(26) expresses the usual relation between the pressure and the value of the pair correlation 
function at contact in hard spheres systems [11]. Note that the low pressure theory therefore 
predicts a (troncated) Gaussian distribution of the overlaps between the spheres while the 
sign of r appearing in (25) may become negative at sufficiently high pressures. This happens 
when many spheres are touching themselves. The probability V of their mutual overlap 
being equal to Q (or —Q) it then becomes higher than the probability V(x = 0) of finding 
"orthogonal" spheres. 

At low Q, the pressure of this solution diverges as 

P = q^q~ (27) 

^% ^ mm 



where A = and Q mm = J I - -. 



III. NUMERICAL SIMULATIONS - THE GLASS TRANSITION 

In order to locate a possible transition, we have performed simulations with Monte Carlo 
dynamics. The simulation starts from a configuration without overlaps between particles, 
and a large value of Q (corresponding to small particles or, equivalently, low density). The 
particles are moven randomly, and the changes are accepted if they do not violate the hard- 
sphere condition. 
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In a first set of simulations (corresponding to infinite pressure in the transient), every 
At steps an attempt is made to enlarge all the particles by an amount SQ: Q — >■ Q — 6Q . 
If the change in size does not generate an overlap, then it is accepted. When a given target 
value Q target is reached, the system evolves at fixed Q = Q tar get- Different annealing speeds 
correspond, in this set of simulations, to smaller proposed changes in Q tar get an d larger 
intervals between changes of Q tar get- 

A second set of simulations was done at constant pressure as follows: the particles are 
moven at random as before, never accepting overlaps between them. After every motion 
an attempt is made to change the particle size. If the change generates an overlap, then 
it is rejected. If it does not generate an overlap, then it is accepted with the Monte Carlo 
probability associated to an energy A^ 2 P Q (cfr. Eq (4)). At not too high densities, we 
have observed that the number of sweeps needed to thermalize roughly scales as A^ 2 . 

We have done simulations for p = 2 and D = 4, 8, 16, 32 and 64. For dimensions 
up to 16, the system goes at high densities to a regular structure (see below). We have 
done longer simulations for D = 32 and 64. Fig. 1 shows for D = 64 the plot of average 
interparticle 'squared distance' J2 a ^b qlb/^iN — 1) versus particle size Q obtained with the 
infinite-pressure algorithm for four different annealing times (each one half of the previous, 
the longest run consisting of about 10 7 sweeps), together with the analytical curve. One can 
see the departure of the numerical points from the low-density solution at a value of Q < 1. 
The dependence on the annealing velocity becomes more important at higher densities. 
However, we have checked that the correction due to the finite annealing time vanishes as 
the annealing time to a power close to —0.5. The limit of infinitely slow annealing yields a 
value that is different at high densities from the analytic computation for the liquid phase. 

More detailed information is obtained by plotting the histogram of normalized relative 
distances q a b/Q for different particle sizes Q (Fig. 2). As the particle size is increased the 
tipical distances become more and more concentrated near the minimal distance Q, i.e. 
most spheres are almost touching. 

The low-density solution of Section II predicts a Gaussian form for this distribution in 
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FIG. 2. Histogram of q ab vs. q/Q for Q of 0.791960, 0.862670, 0.919239 and 1.06066 
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0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 

FIG. 3. f vs Q (full line) analytical result for the low-density phase. Coefficient in the Gaussian 
fit of the histogram of q ab from the simulation. In the inset the error corresponding to the Gaussian 
fit. 

the low-density phase. By comparing this prediction with the numerical data we can see at 
which density the low-density solution breaks down. 

In figure 3 we plot the coefficient of the Gaussian obtained by doing a numerical fit of 
the data versus the analytical prediction of the previous section. The agreement is good up 
to a value of Q around 0.84. A further confirmation of the transition density is obtained by 
computing the error in the Gaussian fit. The inset of Fig. 3 shows that the error is negligible 
for Q > 0.84 and starts growing for smaller values of Q (corresponding to the high-density 
phase). This means that, for Q < 0.84, the distribution P(q ab ) is definitely not a Gaussian. 

With the constant pressure algorithm we have performed a slow annealing in pressure. 
The result of Q vs. P for D = 32, p = 2 (the volume versus pressure curve for this model) is 
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5 10 15 20 25 

P 

FIG. 4. Sphere diameter Q as a function of the pressure p. Diamonds : pressure estimated 
from the histogram P(q ab ) for D = 64 using the longest run. Crosses : annealing in the pressure 
for D = 32. Full line : low-density analytical calculation. 

shown in Fig. 4. Again, the analytical low-pressure curve breaks away from the numerical 
one at a value of Q consistent with the previous one. In that figure we also show the pressure 
as obtained from the probabilty of q a b at q a b = Q (see Ref. [11]). An asymptotic form of the 
pressure (27) describes well the data for the longest run at low Q and D = 64, the values 
of A and Q m in being 0.4 and 0.76 respectively. These values are close but different from the 
theoretical low density predictions of A = * and Q m in = 

From these simulations we see that there is some sort of dynamical freezing at large 
pressures. In order to try to undestand the nature of the dynamical high-pressure regime, 
we have performed a numerical 'aging experiment' at P = 8.4 (D = 32, p = 2). We perform 
a rapid 'quench' in pressure down to P = 8.4 in the presence of a constant field (h = 30) * 
as in equation (5). The conjugate density Eq. (6) tends to stabilize after a transient. We 
now cut the field at different times t w , and observe the decay of the density response (6) for 
further times t. The results are shown in Fig. 5, for three different t w , together with (in the 
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FIG. 5. Density response in the aging simulation m(t + t w ) vs. /, log-linear plot. The waiting 
times are t w = 25. 10 4 , 63. 10 4 and 102. 10 4 . Inset : density response for constant field; the vertical 
range is from 0.15 up to 0.20. The time axis is the same for both graphs. 

inset) the evolution of the density reponse for constant field. 

We see that the remanent density response has a very slow decay after the field has been 
cut off. The spheres are very blocked in a position that was generated under the field, and 
have difficulty in rearranging themselves. The system in which the field has been cut off 
first has the fastest decay of the density response. The question now is whether the system 
has 'true' or 'weak' ergodicity breaking, as we shall discuss below. 

IV. DISCUSSION 

Let us now discuss the physics of this model, in view of the results we have obtained and 
of what we know from spin-glass dynamics. First of all we must discuss possible finite-size 
effects. A finite system cannot have a sharp transition; at finite pressure it will eventually 
equilibrate, and thus cannot have a glassy regime either. We have simulated the system with 
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D = 16, p = 2 and found evidence of equilibration: the distribution of overlaps tends to have 
a few peaks, suggesting a fall into a 'crystalline' state. Though we have not performed a 
detailed analysis of the dependence of the equilibriation time upon D and iV, we know 
that it rapidly increases with both parameters. For the larger sizes we have presented {e.g. 
D = 32, 64, p = 2) we have observed no evidence for the formation of a crystalline state, 
the resulting function P(q ab ) being smooth. 

Let us then discuss what is the behaviour of the model for large but finite times (as 
compared to the system size). The first question that arises is what is the role of the 
annealing procedure or, rather, what would it be in a system with N — > oo. In order 
to distinguish the dependence upon different annealing velocities from aging effects, we 
shall consider first the asymptotic value of one-time quantities (distribution of interparticle 
distances, density, etc). Let us be more precise and consider the following procedure: we 
start with the system at low-density and we increase the pressure with a given velocity up 
to a given final pressure. After that we let the system relax at constant pressure for a long 
subsequent time. The question is now whether the asymptotic value of, say, the distribution 
of interparticle distances depends upon the annealing velocity. Our data are compatible 
with the assumption that there is no such dependence. In this respect this model behaves 
as real spin-glasses which are different from ordinary glasses in that the large time limit of 
one-time quantities are independent of the annealing procedure [15]. 

In the high-pressure regime there is obviously some form of slow dynamics (see Fig. 
5). The question now arises as to whether we are in the presence of 'strong' or (as is 
the case in spin glasses) 'weak' ergodicity breaking. In the first case, the system remains 
confined in a small region of phase-space and remembers forever the conditions in which it 
was prepared. In the 'weak' ergodicity breaking scenario, there is no confinement within a 
small region of phase space. For generic initial conditions after a sufficiently large time the 
system forgets any perturbation that lasted for a finite time (and in particular the initial 
condition). However, as the system ages it moves slower and slower, and the time needed 
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to forget any perturbation grows with the time during which the perturbation acted. 

In the context of our 'aging experiment' of the Section III, the questions are whether the 
density response decays to zero for any t w (in which case the ergodicity-breaking is 'weak') 
and whether the typical decay time depends on t w ('aging'). 

The results shown in Fig. 5 seem to point in the direction of weak ergodicity with aging, 
but they are hardly conclusive. The decay is extremely slow, an effect that has already been 
observed in the models with self-induced disorder that appeared recently in the litterature. 
For the related matrix model of Ref. [8] it is clear that the the auto-correlation function 
decays to zero but it does so extremely slowly. For t w 3> 1 and t w /r + t w = 0(1) : 

C(T + t W} t w ) ~ (t w /T + t w )°-° 5 . 

In this note, we have presented a spherical mean-field model for hard-spheres whose 
dynamics is potentially solvable in the high pressure phase. We believe that it can provide 
an analytical description of such typical glassy phenomena as the viscosity increase around 
the glassy transition and aging effects [16] but possibly not remanent dependences on the 
annealing procedures. 
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